---
title: "Saliency_EEG"
output: html_document
date: "2024-09-19"
---

```{r, include=FALSE}
library(base)       # The R Base Package
library(rstudioapi) #to set directory to current file path
library(matrixStats)# Functions that Apply to Rows and Columns of Matrices (and to Vectors) 
library(sjmisc)     # Data and Variable Transformation Functions 
library(stringr)    # Simple, Consistent Wrappers for Common String Operations (filtering of partial name)
library(purrr)      # Functional Programming Tools

library(ggplot2)    # Create Elegant Data Visualizations Using the Grammar of Graphics
library(ggpubr)     # 'ggplot2' Based Publication Ready Plots
library(ggsignif)   # Significance Brackets for 'ggplot2'
library(ggsci)      # Scientific Journal and Sci-Fi Themed Color Palettes for 'ggplot2'
library(ggrain)     # A Rainclouds Geom for 'ggplot2'
library(plotly)     # smoothens the lines in plots, especially the horizontal ones
library(RColorBrewer) # color palette 
cbPalette <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")# color-blind friendly palette

library(psych)      # Procedures for Psychological, Psychometric, and Personality Research
library(irr)        # Various Coefficients of Interrater Reliability and Agreement
library(rstatix)    # Pipe-Friendly Framework for Basic Statistical Tests 
library(broom)      # Convert Statistical Objects into Tidy Tibbles 
library(stats)      # The R Stats Package
library(lme4)       # Linear Mixed-Effects Models using 'Eigen' and S4 
library(lmerTest)   #get p-values for LMM's
library(emmeans)    # Estimated Marginal Means, aka Least-Squares Means
library(car)        # Companion to Applied Regression
library(pbkrtest)   # Parametric Bootstrap, Kenward-Roger and Satterthwaite Based Methods for Test in Mixed Models
library(effectsize) # Indices of Effect Size

library(sjlabelled) # to generate LMM output table 
library(sjPlot)     # table functions
```

### set working directory
```{r}
current_path = rstudioapi::getActiveDocumentContext()$path 
setwd(dirname(current_path ))
print( getwd() )
```

### load data
```{r, include=FALSE}
theta <-read.delim("theta saliency.csv", header=TRUE, sep=";")
alpha <-read.delim("alpha saliency.csv", header=TRUE, sep=";")
beta <-read.delim("beta saliency.csv", header=TRUE, sep=";")
phaseLocked <-read.delim("phase-locked saliency updatedCluster2.csv", header=TRUE, sep=";")

#phaseLocked <-read.delim("phase-locked aggregated.csv", header=TRUE, sep=";") #for exploratory analysis
#alpha <-read.delim("alpha aggregated.csv", header=TRUE, sep=";") #for exploratory analysis

```

### sort table

```{r}
theta_full <-convert_as_factor(theta, vars=  c("latency", "condition", "electrode","subject","frequency_band"))
alpha_full <-convert_as_factor(alpha, vars=  c("latency", "condition", "electrode","subject","frequency_band"))
beta_full <-convert_as_factor(beta, vars=  c("latency", "condition", "electrode","subject","frequency_band"))
phaseL_full <-convert_as_factor(phaseLocked, vars=  c("latency", "condition", "electrode","subject","frequency_band"))
```


## check assumptions / outliers

```{r}
shapiro.test(phaseL_full$amplitude)

check_data <- phaseL_full
lmm <- lmer(amplitude ~ condition*latency+ (1|subject), data=check_data)
qqnorm(resid(lmm))
qqline(resid(lmm))

check_data$cook <- cooks.distance(lmm)
plot(check_data$cook, check_data$subject)

#log transformation to make data more conform with normality (if necessary)
check_data$amp_log <- log(check_data$amplitude*1+1)
lmm_log <- lmer(amp_log ~ condition*latency+ (1|subject), data=check_data)
qqnorm(resid(lmm_log))
qqline(resid(lmm_log))

check_data$cook <- cooks.distance(lmm_log)
plot(check_data$cook, check_data$subject)

check_data$amplitude <- replace(check_data$amplitude, c(check_data$cook > 1),NA) #replaces amplitude with NA for all rows in which cook's distance is larger than 1
check_data$amp_log <- replace(check_data$amp_log, c(check_data$cook > 1),NA)

phaseL_full <- check_data

```

```{r}
ggboxplot(phaseL_full, x="latency", y="amplitude", facet.by = "condition",add=c("mean_se"), title= "phase-locked", xlab= FALSE, ylab="Amplitude [uV]")  + geom_exec(geom_line, phaseL_full, x="latency", y="amplitude", group="subject", color="grey" )+ stat_compare_means(method="t.test")

# test difference between baseline and oddball for each cluster /condition
phaseL_full %>%
  group_by(electrode) %>%
  wilcox_test(data =., amplitude ~ latency, paired=TRUE) %>%
  adjust_pvalue(method = "bonferroni") %>%
  add_significance("p.adj")

phaseL_full %>%
  group_by(electrode) %>%
  t_test(data =., amp_log ~ latency, paired=TRUE) %>%
  adjust_pvalue(method = "bonferroni") %>%
  add_significance("p.adj")

#calculate relative change at oddball frequency relative to baseline
PL_highO <- filter(phaseL_full, condition=="high oddball")
PL_lowO_cl1<-filter(phaseL_full, condition=="low oddball" & electrode=="clusterCPz-Cz-C2-FCz-FC2-Fz")
PL_lowO_cl2<- filter(phaseL_full, condition=="low oddball" & electrode=="clusterC3-C5-CP5") 

diff_ampPL_HO <-as.data.frame( (filter(PL_highO, latency=="aggregatedH 0.125"))$amplitude - (filter(PL_highO, latency=="0.5"))$amplitude)
colnames(diff_ampPL_HO) <- "rel amp HO"
diff_ampPL_LO1 <-as.data.frame( (filter(PL_lowO_cl1, latency=="aggregatedH 0.125"))$amplitude - (filter(PL_lowO_cl1, latency=="0.5"))$amplitude)
colnames(diff_ampPL_LO1) <- "rel amp LO1"
diff_ampPL_LO2 <- as.data.frame((filter(PL_lowO_cl2, latency=="aggregatedH 0.125"))$amplitude - (filter(PL_lowO_cl2, latency=="0.5"))$amplitude)
colnames(diff_ampPL_LO2) <- "rel amp LO2"

rel_amp_all <- stack(cbind(diff_ampPL_HO,diff_ampPL_LO1,diff_ampPL_LO2))

#compare relative change between conditions / clusters
rel_amp_all %>%
  wilcox_test(data =., values ~ ind, paired=TRUE) %>%
  adjust_pvalue(method = "bonferroni") %>%
  add_significance("p.adj")

relAmpPLot <- ggerrorplot(rel_amp_all, x="ind", y="values", desc.stat= "mean_ci", combine = TRUE, size= 1, 
          add=c( "jitter"),add.params = list(color="black", alpha=0.5), title= FALSE, xlab= FALSE, ylab="∆ amplitude [µV]", color="ind", palette = c("Red", "Blue","Blue"), legend="RIGHT", panel.labs.background=list(color="white", fill="white"), panel.labs.font=list(face="bold", size=13))+
  theme(panel.border = element_blank(), axis.line = element_line()) + 
    scale_x_discrete(labels=c('high oddball', 'low oddball cluster 1', 'low oddball cluster 2'),guide = guide_axis(n.dodge = 2))+
    stat_compare_means(method="wilcox.test", paired=TRUE, comparisons = list(c("rel amp HO", "rel amp LO1"),c("rel amp HO", "rel amp LO2"))  ,label = "p.signif", tip.length = 0.02)#+
    stat_summary(data= peaks_highO_rel,  fun = "mean", geom = "line", aes(group = subject),color="grey", alpha=0.5)

```

```{r}
ggerrorplot(rel_amp_all, x="ind", y="values", desc.stat= "mean_ci", combine = TRUE, size= 1, 
          add=c( "jitter"),add.params = list(color="black", alpha=0.5), title= FALSE, xlab= FALSE, ylab="∆ amplitude [µV]", color="ind", palette = c("Red", "Blue","Blue"), legend="RIGHT", panel.labs.background=list(color="white", fill="white"), panel.labs.font=list(face="bold", size=13))+
  theme(panel.border = element_blank(), axis.line = element_line()) + 
    scale_x_discrete(labels=c('high oddball', 'low oddball cluster 1', 'low oddball cluster 2'),guide = guide_axis(n.dodge = 2) ,expand=expansion(mult=c(0,1)))+
    stat_compare_means(method="wilcox.test", paired=TRUE, comparisons = list(c("rel amp HO", "rel amp LO1"),c("rel amp HO", "rel amp LO2"))  ,label = "p.signif", tip.length = 0.02)#+
    stat_summary(data= peaks_highO_rel,  fun = "mean", geom = "line", aes(group = subject),color="grey", alpha=0.5)
    
ggarrange(relRatPlot, relAmpPLot,  common.legend=TRUE,legend= "none", labels = c("A","B"), ncol = 2, nrow=1,align = "hv") #reRatPlot has to be imported from "saliencyVAS.Rmd" file, illustration of relative difference in ratings across conditions
```
